NPSEC-93-020 



NAVAL POSTGRADUATE SCHOOL 

Monterey, California 




Analysis Using Bi-Spectral 
Related Technique 

by 

Ralph Hippenstiel 

// 

November 17, 1993 



Approved for public reslease; distribution unlimited. 

ired for: Naval Command Control & Ocean Surveillance Center 

Code 732 
271 Catalina Blvd. 

San Diego, CA 92152-5000 



FedDocs 
D 208.14/2 
NPS-EC-93-020 



O >0 ! 

,?J/. !* 






Naval Postgraduate School 

Monterey, California 93943-5000 



Rear Admiral T. A. Mercer 
Superintendent 

This report was funded by the Naval Command Control and Ocean Surveillance Center. 
Reproduction of all or part of this report is authorized. 

This report was prepared by: 



H. Shull 
Provost 



REPORT DOCUMENTATION PAGE 


Form Approved 
OM8 No 0704-0188 


PuDlic reporting Durden tor this collection of nformanon is estimated to average ' hour oer response including the time ♦or reviewing instructions, searcrtmg e*istmg data sources, 
gathering and maintaining the data needed, and comoietmg and reviewing tne collection of information Send comments regarding this Durden estimate or any other avoect of this 
collection of information, including suggestions tor reducing this Duroen to Washington neadauarters Services. Directorate for information Operations and Reports 1 2 1 5 Jefferson 
Davis Highway. Suite 1204 Arlington, j A 22202-4302 and to the Office of Management and Sudget. Paperwork Reduction Protect (0704-0 T 88 ). Washington. DC 20S03 


1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

November 17, 1993 Final Report Oct 1. 92-Sep 30, 93 


4. TITLE AND SUBTITLE 

Analysis Using Bi-Spectral Related Techniques 


5. FUNDING NUMBERS 


6. AUTHOR(S) 

Ralph Hippenstiel 


7. PERFORMING ORGANIZATION NAME(S) AND AODRESS(ES) 

Naval Postgraduate School 
Monterey, CA 93943-5000 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

NPSEC-93-020 


9. SPONSORING /MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

Naval Command Control L Ocean Surveillance Center 

Code 732 

271 Catalina Blvd. 

San Diego. CA 92152-5000 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


11. SUPPLEMENTARY NOTES 

The views expressed in this report are those of the authors and do not reflect the 
official policy or position of the Department of Defense or the United States Government. 


12a. DISTRIBUTION /AVAILABILITY STATEMENT 

Approved lor public release; 
distribution unlimited. 


12b. DISTRIBUTION CODE 


13. ABSTRACT (Maximum 200 words) 



This report investigates the use of Bi-spectral related techniques to extract detection/classification clues from time- 
frequency representations. Earlier results have indicated that 1.5-D spectral techniques (a degenerate version of the 
Bi-spectrum) has potential Signal to Noise Ratio (SNR) gain over conventional techniques. This is partially due to 
the rejection of Gaussian like perturbations by the cumulant based techniques. The third order moment for Gaussian 
zero mean random processes is essentially zero. It is assumed that i) the independence of the Instantaneous Power 
Spectrum (IPS) from analytic signal representations and ii) the Gaussian character of the noise causes minimal 
distortion of the output variables, hence permit signal characterization at low to moderate SNRs. The performance 
is demonstrated on i) synthetic signals and ii) on segment(s) of the NOSC data set. 



14. SUBJECT TERMS 

spectral estimation, instantaneous spectrum. Bi-spectrum, 
cumulant based spectrum, detection 


15. NUMBER OF PAGES 

60 


16. PRICE CODE 


17. SECURITY CLASSIFICATION 
OF REPORT 

UNCLASSIFIED 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

UNCLASSIFIED 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

UNCLASSIFIED 


20. LIMITATION OF ABSTRACT 

SAR 



\$N 7540*01-280-5500 . Standard Form 298 (Rev 2-89) 

i Prescribed Dy ANSI Std 239- ’8 



Analysis Using Bi-Spectral Related Techniques. 



Ralph Hippenstiel 

Electrical and Computer Engineering Department 

Naval Postgraduate School 
Monterey, CA 93943-5000 



November 17, 1993 

Final Report 
on Research sponsored by 
Naval Ocean System Center 
San Diego, CA 



1 



Analysis Using Bi-Spectral Related Techniques. 



Table of Content. 



i. Objectives. 

1. Introduction. 

2. Instantaneous power spectrum (IPS), cumulant based 
instantaneous power spectrum (CUM) , and spectrogram based 
representation (LOFAR) . 

2.1 IPS 

2.2 CUM 

2.3 LOFAR 

3. Theoretical results for averaged surfaces. 

3 . 1 Averaged IPS 

3 . 2 Averaged CUM 

3 . 3 Averaged LOFAR 

3.4 Deflection ratio comparison 

4. Simulation results using sinusoids embedded in white Gaussian 
noise . 

5. Processing results using NOSC data base. 

6. Conclusion and recommendation. 

7. Program listings. 

8. References. 

9. Appendix A. 



2 



i. Objectives: 



This report investigates the use of Bi-spectral related 
techniques to extract detection/classification clues from time- 
frequency representations. 

Earlier results have indicated that 1.5-D spectral techniques (a 
degenerate version of the Bi-spectrum) has potential Signal to 
Noise Ratio (SNR) gain over conventional techniques. This is 
partially due to the rejection of Gaussian like perturbations by 
the cumulant based techniques. The third order moment for Gaussian 
zero mean random processes is essentially zero. It is assumed that 
i) the independence of IPS from analytic signal representations 
and ii) the Gaussian character of the noise causes minimal 
distortion of the output variables, hence permit signal 
characterization at low to moderate SNR's. 

The performance is demonstrated on i) synthetic signals and ii) on 
segment (s) of the NOSC data set. 



1. Introduction: 

This report examines the instantaneous power spectrum (IPS) 
and the cumulant based version of IPS (CUM) to address the question 
whether or not these techniques offer possible improvements 
relative to conventional spectral processing (i.e. spectrogram, 
Lofargram, periodogram) . 

Section 2 provides the mathematical definition for IPS, CUM, 
and LOFAR while section 3 allows a theoretical comparison of these 
techniques. The comparison is performed using a deflection 
criterion (i.e. maximum signal magnitude versus standard deviation 
of the detector output under noise only conditions) . 

Section 4 contains simulation results, while section 5 
provides actual SONAR data results. 

Section 6, and 7 contain the conclusion and program listings, 
respectively. Section 8 provides the list of references. 
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2. IPS, CUM and LOFAR 

2.1 IPS: 

IPS is a member of the Cohen's class [1], and is also known as 
the real part of the Rihaczec distribution. Our particular 
implementation differs in that filtering is employed as the data is 
processed (equation 1) . Earlier results have shown that in contrast 
to the Wigner-Ville Distribution (WVD) no spectral cross modulation 
terms are generated. However, it has wider spectral peaks and 
superimposed spectral auto modulation terms [2]. These auto terms 
do not interfere with the interpretation of T-F surfaces, but can 
help in the detection of a weak stable spectral line component. The 
digital implementation of IPS is given by 

IPS{n,k) (x{n) x" (n-m) + x* (n) x(n+m) ) w(m) exp-j2nkm/N 

eq. 1 

where k = frequency index, n = time index, N = the length of data 
used, m = shift parameter and w( ) = lag window. 

Some earlier results are published in [2,3,4]. 

2.2 CUM: 

In general the Bi-spectrum is defined as 

(ii 2 ) ='ZiXh Cx{1 ' ,l2) 

where the cumulant is given by [5] 

C x (l 1 ,l 2 ) = E{X' [n)X{n+l r )X{n+l 2 ) ) . 
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A degenerate version of the cumulant is given by 



C X (0,1 2 ) = E(X' (n)X(n)X{n + l 2 ) ) = E( \X(n) \ 2 X(m 1 2 ) ) . 



When replacing the expected value with the instantaneous value 
the degenerate Bi-spectrum, also called the lh D spectrum, becomes 

S x (co) = lx(n) l 2 x(n + l) e~ j<J>1 . 

The cumulant version of IPS is given by 
CUM(n,k ) = — f2 ~\ {\x{n) \ 2 x" {n-m) + \x(n) \ 2 x(n+m) ) w(rn) exp -j 2n km/ N 

2 -N/ 2 

eq. 2 

Cross terms at locations where the true spectrum has no support 
are not created. Also contrary to the conventional Bi-spectrum none 
of the spectrally related components are suppressed. 

In part of the work in [6], the approach of Petropulu [7] which 
uses an instantaneous higher order moment slice was utilized. This 
approach works well, but does not permit separation of spectral 
components having the same spectral dynamics. That is, all 
stationary spectral components (i.e. stable line components) map to 
the same detection cell, making this technique useless in 
scenarios with several similarly behaving spectral components. 

2.3 LOFAR: 

The LOFAR algorithm uses a Hamming windowed Fast Fourier 
Transform (FFT) and displays the magnitude of the transform. 
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To keep T-F dimensions relative small, no zero padding is used. 



P x (n, k) 



E n+N 
m=n 



+AT-1 



x{m) e' j2TlkTn/N | 



eq. 3 

The IPS, CUM, and LOFAR surfaces are averaged (in magnitude) over 
some segments of the time axis leaving only a display of spectral 
magnitude versus frequency. For each segment one spectral line is 
created resulting again in a time frequency display. The averaging 
technique is also known as power or incoherent averaging. 

3. Theoretical results for averaged surfaces. 

The analysis in this section is done by disregarding any 
window effects and assuming that under the noise only hypothesis 
(H 0 ) Gaussian white noise and under the signal only hypothesis (H^ 
a deterministic signal of the form x(n)=A cos ( 27rk 0 n/N) is present. 
The value of k Q , n, and N are assumed to be integers which ensures 
that the spectrum peaks at a bin center of an N-point Fourier 
transform. In all detection schemes (averaged IPS, averaged CUM, 
and averaged LOFAR) magnitudes are averaged since they are 
numerically more conveniently computed then the magnitude squared 
values. 

The comparison is done for normalized versions of 
I (k) = E n | IPS (n, k) | , 

C (k) = S n | CUM (n, k) | , and 

P(k) = Z n |windowed Fourier transform (x(n 1 ,k))| ; where n 1 is the 

time counted from reference point n. 

In what follows the time (n) and frequency (k) dependency is 
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usually suppressed. 



3.1 Averaged IPS. 



— -> 


real 




1 1 
1 1 


r— > 


1/M S n 






y z w 

where x(n) ~ N(0,u 2 ) i.i.d. and M is the appropriate number 
degrees of freedom. 



of 



Under H 0 : 

Under noise only condition (H Q ) the random variable Y denoted by Y 0 
has the following moments [6]: 

E{Y 0 )= a 2 , E { Y Q ) 2 = [ (N+5) / 2 ] a 4 , var(Y 0 ) = [(N+3)/2] a 4 . 

For the transformation Z = SQRT(Y 2 ), 

the resulting probability density function is given by 
f z (z) = (f Y (z) + f Y (-z) ) U (z) . 

Making the assumption that for large enough number of data points 
the random variable Y is approximately Gaussian distributed we can 
derive the moments for Z as 
E(Z„) = SQRT{ (N+3)/167T} a 2 e' (2/(N+3)) 

E(Z 0 2 ) = E { Y Q ) 2 = ( (N+5 ) /2 ) a 4 
var(Z 0 ) = [(N+5J/2 -({(N+3) e ‘ (4/(N+3)) }/ 16 t r) ] a 4 
var(W Q ) a 1/M var(Z 0 ) (for large N) 
a (1/M) 0.48 N a 4 . 

A typical surface has about N/3 degrees of freedom [6], providing 
a standard deviation (under the H 0 hypothesis) of 
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o(W 0 ) ~ 1.2 a 2 . 

Under H, (at spectral location ± k Q ) : 
Z 1 = |yJ = (NA 2 ) cos 2 (27rk 0 n/N) 

W, = (NA 2 )/M S n cos 2 (27rk 0 n/N) 

« (NA 2 )/2. 



The deflection D ps , for IPS becomes 

D psi = W 1 /a(W 0 ) = (NA 2 /2 ) / ( 1 . 2 a 2 ) = (N/2 . 4 ) ( A 2 /a 2 ) . 



3.2 Averaged CUM. 



CUM 




real 




1 1 
1 1 




l/M E„ 




c 




Y 


Z 





Under H 0 : 

E(C 0 ) = 0 + jO 

E(C 0 ) 2 = var(C 0 ) 

= (6N+54)/4 (a 2 ) 3 = var(Y 0 ) ; [6] 

(since real part of cum is zero) . 

Z 0 = I Y o I 

f 2 (z)= 2/ (27ra Y 2 ) 1/2 exp- ( z 2 / (2a y 2 ) ) U(z) 

E(Z 0 ) = SQRT ( 2/7T ) a y 

E ( Z Q ) 2 = 3/2 a Y 2 

var(Z 0 ) = (3/2 - 2A) a Y 2 

E(W 0 ) = E(Z 0 ) = SQRT ( 2/tt ) a y 

E(W q ) 2 = 1/M { E ( Z Q ) 2 + (M-l) E 2 (Z q ) ) 
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var(W 0 ) = 1/M { (3/2-2/7T) (6N+54)/4 } a 6 
a(W Q ) = SQRT { 1/M ( 3/2-2 /jt) ( 6N+54 ) }/2 a 3 . 

A typical surface has about N/3 degrees of freedom [6] , so the 

standard deviation becomes 

<7(W 0 ) ss SQRT { 3/N ( 3/2-2/n) ( 6N+54 ) }/2 a 3 . 

Under H t (at spectral location ± k Q ) : 
x(n) = A cos (27rk 0 n/N) 

Y 1 = NA 3 /2 cos 3 (27rkn/N) @ k = ± k Q 

Z 1 = NA 3 /2 | cos 3 ( 27rkn/N) | @ k = ± k Q 

From numerical evaluations we know that the summation is lower 
bounded by 

— V'" -1 I (cos 3 (2nln/M) ) I > 0.424 . 

0 1 ' 

Hence Z 1 > 0.424 NA 3 /2. 

The deflection D C(JM becomes 

D cuh = (0*424 NA 3 / 2 ) SQRT (M) /{ SQRT [( 3/2-2/ir )( 6N+54 ) /4 ] a 3 } 

~ 0.1075 N (A /a) 3 . 



3.3 LOFAR: 



x(n) 




X(k) 



Under H 0 : 

X (k) = x r (k) + j X ; (k) 
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n+N-1 



x(m) exp- j ( 27rkm/N) ; where n is the reference time 



m=n 

E(X r (k))= E(x r )= E (Xj (k) ) = E (Xj ) = 0 
E(X r ) 2 = var(X p ) = var(Xj) = o 2 N/2 
Z = SQRT(X r 2 + X^ ) 

E(Z 0 ) =SQRT(7r/2) SQRT (N/2 ) a 

var(Z 0 ) = ( 2-7T/2 ) a 2 N/2 

E(W 0 ) = E(Z 0 ) =SQRT(7t/ 2) SQRT (N/2) a 

E(W 0 ) 2 =1/M { E ( Z Q ) 2 + (M-l) E 2 ( Z Q ) ) 

var(W Q ) = 1/M var(Z Q ) = (1/M) (2 -tt/ 2 ) a 2 (N/2 ) 

a(W Q ) = SQRT { ( 1/M) (2-7T/2) (N/2)) a 

Under H, (at spectral location ± k Q ) : 

X, (k) =A N/2 @ k = ± k 0 

Z, = A N/2 @ k = ± k 0 

E (W, ) = E ( Z, ) = A N/2 
The deflection D lf .,. D becomes 

LOFAR 

d lofar = (AN/2) /{ (1/M) (2-7T/2)a 2 (N/2) } 1/2 
= N 1/2 [3 A 2 / ( ( 4 — 7T ) cr 2 ) ] 1/2 . 



3.4 Deflection ratio comparison: 

We can now plot the deflection ratio, that is compare and 
determine which detector is more advantageous under which 
conditions . 

Comparing D, ps and D L0FAR : 

D IPS > d lofar 
simplifies to 

N > 20.13 (a 2 / A 2 ) 

Comparing D CUH and D L0FAR : 
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^CUM > ^LOFAR 

simplifies to 
N > 302.4 (a 4 / A 4 ) . 

We notice that the processing length variable N, for IPS relative 
to LOFAR has a (SNR)' 1 and for CUM relative to LOFAR has a (SNR)' 2 
dependency . 

4. Simulation results for sinusoid in white Gaussian noise. 

Ten test sets, of length 4096, consisting of sinusoidal 
signals embedded in white Gaussian noise are processed by IPS, CUM 
and LOFAR. Each test set contains seven (7) sinusoids whose SNR 
goes from -3 dB to -21 dB in 3 dB steps. The decrease is monotonic 
relative to the spectral locations so that the SNR decreases in 3 
dB steps as the frequency increases. The sinusoidal frequencies 
are fixed in all test sets at digital frequencies: 

{G>} = { 1 1 Vs , 171/3, 2 9 Va , 37%, 59%, 77%, 111% } * 2?r/2 56 . The phases 

within each set and over all set are independently distributed 
(i.e. no harmonic relationship between different spectral 
components and different realizations) . We note that none of the 
sinusoids has an integer number of cycles over 4096 data points 
ensuring all signals will have side lobes if processed by an FFT of 
size 4096 or less. To demonstrate potential gain the deflection 
ratio as defined in section 3 is computed for each sinusoid and 
averaged over the 10 realizations. 

Figures 1, 2, 3, 4, 5 and 6 are plots of the outcomes of the 
simulation for 128, 256, 512, 1024, 2048 and 4096 data points, 
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respectively. The top part of each figure superimposes all of the 
10 realizations while the bottom part shows the average of the 10 
realizations. We note for example, that the weakest sinusoid (bin 
223-224 of 512 point data length) clearly shows up in the CUM and 
IPS representation (figure 3. A) but is not detectable in the LOFAR 
representation (figure 3.B). This illustrates that the alternative 
representations (i.e. averaged IPS or averaged CUM) can complement 
the traditional (averaged LOFAR) representation. Figures 7-12 
provide plots of the experimental deflection ratios versus SNR in 
the top half of the figures. The solid, dashed, and dotted lines 
represent the experimental deflection for CUM, IPS, and LOFAR, 
respectively. The bottom halves are plots of the difference of 
deflection ratios (D ]ps - D C(JM ) . Figure 13 attempts to reconciliate 
theoretical bounds from section 3 with the experimental evidence 
from this section. The straight solid lines demarcate the 
theoretical transitions points at which CUM works better than LOFAR 
(upper straight line) and where IPS works better than LOFAR (lower 
straight line) . We mention again that no window action was 
accounted for in the derivation of these lines. Superimposed are 
(zigzag lines) the results taken from figures 7 trough 12. The 
upper zigzag line shows the demarcation line where the experimental 
evidence indicates IPS will do better than CUM (in terms of 
deflection ratios). The lower zigzag line consists of the points, 
as a function of SNR and processing length, where the experimental 
performance of all three techniques (IPS, CUM and LOFAR) tends to 
coincide. Above this line the CUM and IPS technique outperform the 
LOFAR approach with improvements increasing as the input SNR increases. 
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data size 1024, proc length 1024 
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powers of 2 processing length 



5. Processing Results using NOSC Data Base. 

In this section the data set supplied by NOSC is processed. 
Figure 14 is a LOFARGRAM of the data set r2A_14 which is 4,128,769 
data points long. The time origin for all grams is chosen to be at 
the top of each figure. Figure 15. A and 15. B are LOFARGRAMS 
produced at NPS to verify the fidelity of the data transfer from 
NOSC to NPS. To see representative spectral plots we selected the 
most interesting segment, that is from time sample 1,024,001 to 
1,228,800. This corresponds to the region 508 to 408 along the time 
(vertical) axis of figure 15. A. Initial processing using IPS and 
especially CUM showed heavy modulation of the time frequency 
displays. This is caused by the low frequency spikes observable on 
figure 14 and 15. A (i.e. strong DC like pulsations (about 5800 
sample points apart) . In the top part of figure 16 which displays 
the time series, used in the experiments, these modulation spikes 
are clearly seen. To minimize the unwanted amplitude modulation 
effects on the time frequency plots, the data is soft limited by 
scaling all values larger than +90 by dividing them by 4 and by 
scaling all values less than -90 by dividing them by 2. The 
resulting time series is plotted in the bottom part of figure 16. 

The contour plots of the remainder of this section show now 
little amplitude modulation allowing the approximation to gray 
scale via contour plots to work. On a color monitor, with a 
sufficient number of color levels of quantization, the modulation 
may not be bothersome. Figure 17 displays 49 traces (along the y- 
axis) which corresponds to 49 individual estimates of the spectral 
density. Closer examination reveals some details on the IPS plot 
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not so easily seen on the lower LOFARGRAM plot. In particular we 
refer to 3 areas: 

1) tick mark 38-48 (vertical axis) about bin 700 (horizontal 

axis) , 

2) tick mark 1-33 (vertical axis) about bin 850 (horizontal 

axis) , and 

3) tick mark 1-20 (vertical axis) about bin 1300 (horizontal 

axis) . 

A word of caution however when using the simple minded contour 
(MATLAB) plots: peaks of identical height but of different width 
will appear differently. In appendix A the effects of choosing 
different contour levels is illustrated. However for all plots the 
best choice of level was manually selected by examining several 
plots differing in the number of levels selected and then retaining 
the one with the best features. 

Figure 18. A and 18. B are displays of IPS, CUM and LOFAR with 
the lines created 512 time data points apart. Figure 18. A (as well 
as 19. A and 20. A) displays spectral location 51 through 1500 while 
figure 18. B (as well as 19 . B and 20. B) displays spectral location 
351 through 1500. This allows better interpretation of the 
different spectral regions. 

Figures 18, 19 and 20 have the same processing parameters, 

they differ in the number of spectral lines computed. Figure 18 
uses a shift of 512 (i.e. 8:1 overlap) resulting in 390 lines, 

figure 19 uses a shift of 1024 (i.e. 4:1 overlap) resulting in 197 
lines, while figure 20 uses a shift of 4096 (i.e. no overlap) 

resulting in 49 spectral lines. 
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Figure 21 illustrates performance as the processing length is 
increased to 8102 (shift 8192, no overlap) providing 25 spectral 
lines (bin 51 through 3300 are displayed) . 

6. Conclusion and recommendation. 

As illustrated in section 3 and 4 processing gain can be 
obtained for IPS and/or CUM relative to the LOFAR processing 
technique. This is a function of signal and processing parameters. 
It may allow glimpses of targets which otherwise maybe hidden. 
However, the processing gain when achievable is bought at the 
increase in processing cost. IPS and CUM are actually tools 
designed more for transient type signals but as shown can be used 
to obtain gain on stable narrow band signals. Simplified 
processing was not attempted but deserve future examination. The 
analysis in section 3 is done totally disregarding lag and data 
windows which typically are used. Future analysis should address 
the influence of the window function on processing. Extensions to 
extract phase information is not attempted but should be pursued. 
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ipsav(clipped to +/-90 redata x, 1,4096., 32) bin 51:1450 




lofarm(clipped to +/-90 redata, 4096, 4096, 4096) bin 51:1450 




Fig. 17 



ipsav(rf, 1,4096, 32) bin 51:1500, shft 512 




cumav(rf, 1,4096,32) bin 51:1500, shft 512 




lofarm(rf, 4096, 4096, 4096) bin 51:1500, shft 512 




Fig. 18 . A 
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ipsav(rf, 1,4096, 32) bin 351:1500, shft 512 




cumav(rf, 1,4096,32) bin 351:1500, shft 512 




lofarm(rf, 4096, 4096, 4096) bin 351:1500, shft 512 
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ipsav(rC 1 ,4096,32) bin 351:1500, shft 1024 




cumav(rf, 1,4096,32) bin 351:1500, shft 1024 
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ipsav(rf, 1,4096, 32) bin 51:1500, snft 4096 




cumav(rf, 1,4096,32) bin 51:1500, shft 4096 
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ipsav(rf, 1,4096,32) bin 351:1500, shft 4096 




cumav(rf, 1,4096, 32) bin 351:1500, shft 4096 
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7. Program Listings: 

All algorithms are implemented (including their hard copy 
functions) using MATLAB software (The Mathworks Inc.) . The programs 
can be executed on MS-DOS personal computers (i.e. 386 + co- 
processor or 486 based IBM compatible PC's) or on SUN work 
stations. Due to the long data sequence and large T-F surfaces most 
of the work is done on SUN work stations. The hard copies provided 
are obtained by using the contour plotting function of MATLAB. One 
note of caution is appropriate when interpreting the T-F plots. Due 
to the way contour lines are drawn, equal height spectral peaks can 
appear different if they differ in T-F spread, that is a narrow 
peak will be darker than a wide peak even though these peaks are of 
the same height. This is due to the density of contour lines which 
is higher for narrow peaks than it is for wide peaks. Better 
results are obtained when displaying the data on gray scale or 
color scale output devices. The new 4.0 MATLAB version supports the 
function 'image' which on a gray scale or color monitor will 
improve the readability of LOFARGRAMs. 

The color output (using MATLAB 3.5i) on a SUN SPARC-2 work 
station shows superior detail when compared to the hard copies 
supplied in this report. The hard copies are obtained using a 
screen dump routine which preserves less detail than a postscript 
derived output. The postscript output is not used due to the length 
of time involved generating it on the processing system. 
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>> atypical setup for obtaining IPS and CUM like surfaces 
>> %for i=start:finish 

>> %p(i , :)=ipsav(rf ((i-1)*shft+1 : (i-1)*shft+N),1 ,N,step); 

>> %q(i ,:)=cumav(rf ((i-1)*shft+1 : (i-1)*shft+N),1 ,N,step}; 

>> % N= proc length, step=stepsize, N/step=number of terms averaged to 
>> ^create one spectral line, shft = amount of data shifted to create 
>> %the next surface which will be averaged to create one spectral 
>> fcline at location i 
>> %end 
>> 

>> atypical setup for obtaining LOFAR like surfaces 
>> %for i=start:finish 

>> %po(i ,:)=lofarm(rf C(i-1)*shft+1 : (i-1)*shft+N),id,st,tl); 

>> %id=number of data points in transform, st=stepsize, 

>> %tl =transform length 

>> % Note: rf = input sequence, 

>> % Note: in all run executed id=st=tl 
>> % So the overlap was controlled with shft 
>> 

>> 

SN ■ 
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% function [P, f reqindex] =ipsav (data, wintype, winlen, step) ; 

%This function will calculate a single averaged Instantan. Power Spectral (IPS) 
I line. 

IThe IPS surface characteristics are determined by the selection of window 
%type (wintype), window length (winlen) and the distance that the window 
%is moved through the data sequence (step) . 

% 

%The P matrix plots only positive frequencies (in magnitude) . The 
%outputs timeindex and freqindex can be used in plots to interpret the 
%results . 

%The inputs are: 

%data - The input data string row vector 

%wintype: '0' Rectangular Window 

1 % '1' Hamming Window 

Iwinlen - The desired width of the window, normally half of the siglen 
%step - Time step desired, normally '1' or a multiple of '2' 

%See also IPSSURF, IPSLOFAR 

function [P, freqindex] =ipsav (data, wintype, winlen, step) 

[datarows, datacolumns ] =size (data) ; 
if datarows ~=1 
data=data . ' ; 

end 

siglen=length (data) ; 
if wintype==0 

win=ones (winlen-1, 1) ; 
elseif wmtype==l 
win=hamming (winlen-1) ; 

end 

W=[win (winlen/2 :-l : 1) ] ; 

x— [ zeros ( 1 , winlen) data zeros (1, winlen) ] ; 

P=zeros (1, winlen/2) ; 

for n=wmlen+l : step : siglen-fwinlen-step+1 

Xm= [con j (x (n : — 1 : n— (winlen/ 2-1) ) ) . ' x (n:n+ (winlen/2- 1 ) ) . ' ] ; 

Xn= [x (n) ; con j (x (n) ) ] ; 
product= ( (Xm*Xn) . *W) . ' ; 

product= [product 0 con j (product (winlen/2 : — 1 : 2 ) ) ] ; 
ptemp=f ft shift (real ( . 5* f ft (product) ) ) ; 

P=P+abs (ptemp (winlen/2 + 1 :winlen) ) ; 
end 

[prow, pcolumn] =size (P) ; 
freqindex= [0 :pcolumn-l] / 
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% function [P, f reqindex, time index] =cumav (data, wintype, winlen, step) ; 

%This function will calculate the 1.5 D Spectral surface. 

%this surface consists of averaged (shorter) cumulant based surfaces 
%in a power or magnitude averaged sense 

%The 1.5 D surface characteristics are determined by the selection of window 
%type (wintype) , window length (winlen) and the distance that the window 
%is moved through the data sequence (step) . 

%The surface is placed sequentially in the P matrix for display. 

%The P matrix plots only the positive half of the spectral plane. The 
^outputs timeindex and freqindex can be used in plots to interpret the 
Iresults . 

%The inputs are: 

%data - The input data string must be in row vector form 

%wintype: '0' Rectangular Window 

% '1' Hamming Window 

%winlen - The desired width of the window, normally half of the siglen 
%step - Time step desired, normally '1' or a multiple of '2' 

%See also ONESURF, ONELOFAR 



function [P, freqindex] =cumav (data, wintype, winlen, step) 

[datarows, datacolumns] =size (data) ; 

if datarows~=l 

data=data . ' ; 

end 

siglen=length (data) ; 
if wintype==0 

win=ones (winlen-1, 1) ; 
elseif wintype==l 
win=hamming (winlen-1) ; 

end 

W= [win (winlen/2 : — 1 : 1 ) ] ; 

x— [zeros (1, winlen) data zeros (1, winlen) ] ; 

P=zeros (1, winlen/2) ; 

for n=winlen+l : step : siglen+winlen-step+1 
Xn= [abs (x (n) ) ^2 ; abs (x (n) ) ^2 ] ; 

Xm= [con j (x (n :-l :n- (winlen/2-1) ) ) . ' x (n :n+ ( winlen /2-1) ) . ' ] ; 
product= ( (Xm*Xn) . *W) . ' ; 

product= [product 0 con j (product (winlen/2 1:2) ) ] ; 
ptemp=f ftshift (real ( . 5*f ft (product) ) ) ; 

P=P+abs (ptemp (winlen/ 2+1 : winlen) ) ; 

end 

[prow, pcolumn] =size (P) ; 
f reqindex= [ 0 :pcolumn-l ] ; 
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I This program computes the LOFARGRAM based on windowed periodograms and displal ys the resu 
lid=no of points in xform, step = shift, tlen=transf orm length 
I [P ] =lofarm (name, id, step, tlen) 
function [P ] =lo farm (name, id, step, tlen) 
clear pow 

lname=input (' name of input file is '); 

[drow, dcolumn] =size (name) ; 
if drow~=l 

name=name . ' ; 

end 

len=length (name) ; 

lid=input(' No. of non-zero data points in transform '); 

%step=input (' shift (step size) in data points '); 

%tlen=input (' transform length '); 

i=fix ( (len-id) /step) +1; 
win=hamming (id) . ' ; 
for ic=l:i 

ppow=abs (fft (name (1+ (ic-1) *step: id+ (ic-1) *step) . *win, tlen) ) ; 
pow (ic, : ) =ppow (1 : tlen/2) ; 
end 
P=pow; 

clear name drow dcolumn len id step tlen i ic win ppow 
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9 . Appendix A 



This appendix provides 3 figures which illustrate the effects of 
level selection when using a simple (no-gray scale , no color ) 
printer. All three figures use a shift of 4096 (i.e. no overlap of 
data between successive spectral lines). Fig A.l, A. 2 and A. 3 show 
IPS, CUM, and LOFAR respectively. The number of quantization 
levels is different on each figure. The label on each figure is 
self explanatory. 
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